In [1]:
import mne
# from mne.preprocessing import (ICA, create_eog_epochs, create_ecg_epochs, corrmap)
# from mne.decoding import cross_val_multiscore
import matplotlib.pyplot as plt
from matplotlib import cm
import seaborn as sns
from fooof import FOOOFGroup
from fooof.bands import Bands
from fooof.analysis import get_band_peak_fg
from fooof.plts.spectra import plot_spectrum
import numpy as np
/tmp/ipykernel_3386360/16807159.py:10: DeprecationWarning: The `fooof` package is being deprecated and replaced by the `specparam` (spectral parameterization) package. This version of `fooof` (1.1) is fully functional, but will not be further updated. New projects are recommended to update to using `specparam` (see Changelog for details). from fooof import FOOOFGroup
In [2]:
# Load data
root_path = '/home/dspserver/ngHAN/imsp/dataset/archive/data_preprocessed_python/'
In [3]:
# Function to load data from each participant file
import pickle
def read_eeg_signal_from_file(filename):
x = pickle._Unpickler(open(filename, 'rb'))
x.encoding = 'latin1'
p = x.load()
return p
In [4]:
import os
labels = []
data = []
# Load data from first 22 files
for file in os.listdir(root_path)[:22]:
trial = read_eeg_signal_from_file(os.path.join(root_path, file))
labels.append(trial['labels'])
data.append(trial['data'])
# Re-shape arrays into desired shapes
labels = np.array(labels)
labels = labels.flatten()
labels = labels.reshape(880, 4)
data = np.array(data)
data = data.flatten()
data = data.reshape(880, 40, 8064)
Extract EEG data¶
In [5]:
eeg_channels = np.array(["Fp1", "AF3", "F3", "F7", "FC5", "FC1", "C3", "T7", "CP5", "CP1", "P3", "P7", "PO3", "O1", "Oz", "Pz", "Fp2", "AF4", "Fz", "F4", "F8", "FC6", "FC2", "Cz", "C4", "T8", "CP6", "CP2", "P4", "P8", "PO4", "O2"])
peripheral_channels = np.array(["hEOG", "vEOG", "zEMG", "tEMG", "GSR", "Respiration belt", "Plethysmograph", "Temperature"])
In [6]:
eeg_data = np.array([[data[i, j] for j in range(len(eeg_channels))] for i in range(len(data))])
eeg_data = eeg_data.reshape(len(data), len(eeg_channels), -1)
print(eeg_data.shape)
peripheral_data = np.array([[data[i, j] for j in range(32, len(data[0]))] for i in range(len(data))])
peripheral_data = peripheral_data.reshape(len(data), len(peripheral_channels), -1)
print(peripheral_data.shape)
(880, 32, 8064) (880, 8, 8064)
In [7]:
# Making info is important for because later we will use it to indicate the channels
info = mne.create_info(eeg_channels.tolist(), ch_types=32*['eeg'], sfreq=128)
info.set_montage('standard_1020')
print(info)
<Info | 8 non-empty values bads: [] ch_names: Fp1, AF3, F3, F7, FC5, FC1, C3, T7, CP5, CP1, P3, P7, PO3, O1, ... chs: 32 EEG custom_ref_applied: False dig: 35 items (3 Cardinal, 32 EEG) highpass: 0.0 Hz lowpass: 64.0 Hz meas_date: unspecified nchan: 32 projs: [] sfreq: 128.0 Hz >
Topomaps for each band of frequencies¶
In [8]:
# Theta band, first trial
evData_th = mne.EvokedArray(eeg_data[0], info)
times = np.arange(0.05, 0.251, 0.04)
evData_th.filter(4, 8, picks="eeg")
evData_th.plot_topomap(times, ch_type="eeg", average=60, time_unit='s', size=3)
Setting up band-pass filter from 4 - 8 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 4.00 - Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 3.00 Hz) - Upper passband edge: 8.00 Hz - Upper transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 9.00 Hz) - Filter length: 213 samples (1.664 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
Out[8]:
In [9]:
# Alpha band, first trial
evData_al = mne.EvokedArray(eeg_data[0], info)
times = np.arange(0.05, 0.251, 0.04)
evData_al.filter(8, 12, picks="all")
evData_al.plot_topomap(times, ch_type='eeg', average=60, time_unit='s', size=3)
Setting up band-pass filter from 8 - 12 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 8.00 - Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 7.00 Hz) - Upper passband edge: 12.00 Hz - Upper transition bandwidth: 3.00 Hz (-6 dB cutoff frequency: 13.50 Hz) - Filter length: 213 samples (1.664 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
Out[9]:
In [10]:
# Beta band, first trial
evData_bt = mne.EvokedArray(eeg_data[0], info)
times = np.arange(0.05, 0.251, 0.04)
evData_bt.filter(12, 30, picks="all")
evData_bt.plot_topomap(times, ch_type='eeg', average=60, time_unit='s', size=3)
Setting up band-pass filter from 12 - 30 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 12.00 - Lower transition bandwidth: 3.00 Hz (-6 dB cutoff frequency: 10.50 Hz) - Upper passband edge: 30.00 Hz - Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz) - Filter length: 141 samples (1.102 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
Out[10]:
In [11]:
# Gamma band, first trial
evData_gm = mne.EvokedArray(eeg_data[0], info)
times = np.arange(0.05, 0.251, 0.04)
evData_gm.filter(30, 63.9, picks="all")
evData_gm.plot_topomap(times, ch_type='eeg', average=60, time_unit='s', size=3)
Setting up band-pass filter from 30 - 64 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 30.00 - Lower transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 26.25 Hz) - Upper passband edge: 63.90 Hz - Upper transition bandwidth: 0.10 Hz (-6 dB cutoff frequency: 63.95 Hz) - Filter length: 4225 samples (33.008 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
Out[11]:
Topomaps of subjects from HAHV - HALV - LAHV - LALV groups¶
In [12]:
# Getting samples from 4 label groups, same subject
ev_data_hahv = mne.EvokedArray(eeg_data[1], info)
ev_data_halv = mne.EvokedArray(eeg_data[14], info)
ev_data_lahv = mne.EvokedArray(eeg_data[6], info)
ev_data_lalv = mne.EvokedArray(eeg_data[9], info)
In [13]:
# Plot the topographies across different frequency bands
def plot_topo_psd(evData):
plt.tight_layout()
evData.filter(4, 8)
evData.plot_topomap(0.15, ch_type='eeg', average=60, time_unit='s', size=3)
evData.filter(8, 12)
evData.plot_topomap(0.15, ch_type='eeg', average=60, time_unit='s', size=3)
evData.filter(12, 30)
evData.plot_topomap(0.15, ch_type='eeg', average=60, time_unit='s', size=3)
evData.filter(30, 63.9)
evData.plot_topomap(0.15, ch_type='eeg', average=60, time_unit='s', size=3)
In [14]:
plot_topo_psd(ev_data_hahv)
Setting up band-pass filter from 4 - 8 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 4.00 - Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 3.00 Hz) - Upper passband edge: 8.00 Hz - Upper transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 9.00 Hz) - Filter length: 213 samples (1.664 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
<Figure size 640x480 with 0 Axes>
Setting up band-pass filter from 8 - 12 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 8.00 - Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 7.00 Hz) - Upper passband edge: 12.00 Hz - Upper transition bandwidth: 3.00 Hz (-6 dB cutoff frequency: 13.50 Hz) - Filter length: 213 samples (1.664 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
Setting up band-pass filter from 12 - 30 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 12.00 - Lower transition bandwidth: 3.00 Hz (-6 dB cutoff frequency: 10.50 Hz) - Upper passband edge: 30.00 Hz - Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz) - Filter length: 141 samples (1.102 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
Setting up band-pass filter from 30 - 64 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 30.00 - Lower transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 26.25 Hz) - Upper passband edge: 63.90 Hz - Upper transition bandwidth: 0.10 Hz (-6 dB cutoff frequency: 63.95 Hz) - Filter length: 4225 samples (33.008 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
In [15]:
plot_topo_psd(ev_data_halv)
Setting up band-pass filter from 4 - 8 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 4.00 - Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 3.00 Hz) - Upper passband edge: 8.00 Hz - Upper transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 9.00 Hz) - Filter length: 213 samples (1.664 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
<Figure size 640x480 with 0 Axes>
Setting up band-pass filter from 8 - 12 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 8.00 - Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 7.00 Hz) - Upper passband edge: 12.00 Hz - Upper transition bandwidth: 3.00 Hz (-6 dB cutoff frequency: 13.50 Hz) - Filter length: 213 samples (1.664 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
Setting up band-pass filter from 12 - 30 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 12.00 - Lower transition bandwidth: 3.00 Hz (-6 dB cutoff frequency: 10.50 Hz) - Upper passband edge: 30.00 Hz - Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz) - Filter length: 141 samples (1.102 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
Setting up band-pass filter from 30 - 64 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 30.00 - Lower transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 26.25 Hz) - Upper passband edge: 63.90 Hz - Upper transition bandwidth: 0.10 Hz (-6 dB cutoff frequency: 63.95 Hz) - Filter length: 4225 samples (33.008 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
In [16]:
plot_topo_psd(ev_data_lahv)
Setting up band-pass filter from 4 - 8 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 4.00 - Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 3.00 Hz) - Upper passband edge: 8.00 Hz - Upper transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 9.00 Hz) - Filter length: 213 samples (1.664 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
<Figure size 640x480 with 0 Axes>
Setting up band-pass filter from 8 - 12 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 8.00 - Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 7.00 Hz) - Upper passband edge: 12.00 Hz - Upper transition bandwidth: 3.00 Hz (-6 dB cutoff frequency: 13.50 Hz) - Filter length: 213 samples (1.664 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
Setting up band-pass filter from 12 - 30 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 12.00 - Lower transition bandwidth: 3.00 Hz (-6 dB cutoff frequency: 10.50 Hz) - Upper passband edge: 30.00 Hz - Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz) - Filter length: 141 samples (1.102 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
Setting up band-pass filter from 30 - 64 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal bandpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 30.00 - Lower transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 26.25 Hz) - Upper passband edge: 63.90 Hz - Upper transition bandwidth: 0.10 Hz (-6 dB cutoff frequency: 63.95 Hz) - Filter length: 4225 samples (33.008 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.0s
Topomaps of subjects from HAHV - HALV - LAHV - LALV groups by FOOOF¶
In [17]:
# Deal with NaN values when the model cannot detect peaks in any given range
def check_nans(data, nan_policy='zero'):
"""Check an array for nan values, and replace, based on policy."""
# Find where there are nan values in the data
nan_inds = np.where(np.isnan(data))
# Apply desired nan policy to data
if nan_policy == 'zero':
data[nan_inds] = 0
elif nan_policy == 'mean':
data[nan_inds] = np.nanmean(data)
else:
raise ValueError('Nan policy not understood.')
return data
In [18]:
# Plot the topographies across different frequency bands
def plot_psd_fooof(evData):
fg = FOOOFGroup(peak_width_limits=[1, 6], min_peak_height=0.05, peak_threshold=2., max_n_peaks=6, verbose=False)
bands = Bands({'theta': [4, 8],'alpha': [8, 12],'beta': [12, 30],'gamma': [30, 64]})
freq_range = [1, 128]
# Calculate power spectra across the the continuous data by MNE
spectrum = evData.compute_psd(fmin=1, tmin=0, tmax=250, method="welch", n_fft=300, n_overlap=150) # <-- modified here for new mne version
spectra, freqs = spectrum.get_data(return_freqs=True)
fg.fit(freqs, spectra, freq_range)
# Plot the topographies across different frequency bands
fig, axes = plt.subplots(1, 4, figsize=(20, 5))
for ind, (label, band_def) in enumerate(bands):
# Extract the power peaks across channels for the current band
band_power = check_nans(get_band_peak_fg(fg, band_def)[:, 1])
# Create a topomap for the current oscillation band
mne.viz.plot_topomap(band_power, evData.info, cmap=cm.viridis, axes=axes[ind], show=False);
axes[ind].set_title(label + ' power', {'fontsize' : 16})
plt.tight_layout()
plt.show()
In [19]:
def plot_psd_peak(evData):
fg = FOOOFGroup(peak_width_limits=[1, 6], min_peak_height=0.05, peak_threshold=2., max_n_peaks=6, verbose=False)
bands = Bands({'theta': [4, 8],'alpha': [8, 12],'beta': [12, 30],'gamma': [30, 64]})
freq_range = [1, 128]
# Calculate power spectra across the the continuous data by MNE
spectrum = evData.compute_psd(fmin=1, tmin=0, tmax=250, method="welch", n_fft=300, n_overlap=150) # <-- modified here for new mne version
spectra, freqs = spectrum.get_data(return_freqs=True)
fg.fit(freqs, spectra, freq_range)
# Check the largest detected peaks within each band
fig, axes = plt.subplots(1, 4, figsize=(20, 6))
for ind, (label, band_def) in enumerate(bands):
# Get the power values across channels for the current band
band_power = check_nans(get_band_peak_fg(fg, band_def)[:, 1])
# Extracted and plot the power spectrum model with the most band power
fg.get_fooof(np.argmax(band_power)).plot(ax=axes[ind], add_legend=False)
axes[ind].yaxis.set_ticklabels([])
axes[ind].set_title('biggest ' + label + ' peak', {'fontsize' : 16})
plt.tight_layout()
plt.show()
In [20]:
plot_psd_fooof(ev_data_hahv)
Effective window size : 2.344 (s)
In [21]:
plot_psd_peak(ev_data_hahv)
Effective window size : 2.344 (s)
In [22]:
plot_psd_fooof(ev_data_halv)
Effective window size : 2.344 (s)
In [23]:
plot_psd_peak(ev_data_halv)
Effective window size : 2.344 (s)
In [24]:
plot_psd_fooof(ev_data_lahv)
Effective window size : 2.344 (s)
In [25]:
plot_psd_peak(ev_data_lahv)
Effective window size : 2.344 (s)
In [26]:
plot_psd_fooof(ev_data_lalv)
Effective window size : 2.344 (s)
In [27]:
plot_psd_peak(ev_data_lalv)
Effective window size : 2.344 (s)
